%based on MLI2 basic properties from Macosko lab collaboration. Patterns
%are the same name.


close all;
clear all;
fclose all;


path1=['C:\Users\Tom�s\Desktop\Temp\'];
dir_to_search=[path1]
formatpattern = fullfile(dir_to_search, '/*.ibw');
dInfoNS = dir(formatpattern);
for i=1:numel(dInfoNS)
    nameCell{i}=dInfoNS(i).name;
end
nameCellSort=natsortfiles(nameCell);% natural order sort file names.


%%  make a logical vector for each pattern type so I can perform operations on each group of files, regardless of how they were acquired.
for i=1:numel(dInfoNS)
    temp=IBWread([dInfoNS(i).folder '\' nameCellSort{i}]);
    
    if    numel(strfind(temp.WaveNotes,'Spontaneous'))>0  %if the Wavenotes contain the pattern for each logical
        spontLogical(i)=1;
    else
        spontLogical(i)=0;
    end
    
    if  numel(strfind(temp.WaveNotes,'vStep;a1:-10'))>0   %for capacitance and input resistance
        vStepLogical(i)=1;
    else
        vStepLogical(i)=0;
    end
    
    if   numel(strfind(temp.WaveNotes,'vStep;a1:-30'))>0  % for Ih
        IhLogical(i)=1;
    else
        IhLogical(i)=0;
    end
    
    if numel(strfind(temp.WaveNotes,'IStep'))>0  % for IF curve
        IFLogical(i)=1;
    else
        IFLogical(i)=0;
    end
    
end

IFLogical=logical(IFLogical);
IhLogical=logical(IhLogical);
vStepLogical=logical(vStepLogical);
spontLogical=logical(spontLogical);

%% Analyze capacitance and input resistance
if sum(vStepLogical)>0
fileNames=nameCellSort(vStepLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    dataMat(:,i)=temp.y;
end

trace= median(dataMat,2);
% this section was taken from the old script, may be useful
Dx=temp.dx; % 1e5 is sampling rate;
DV=-0.01; %10 mV;
t=(0:(numel(trace)-1))*Dx;
figure;
plot(t,trace);
xlabel('Time (s)');
ylabel('I (pA)');
hold on;
plot(0.1,trace(1e4),'r.','markersize',10);
%intStart= find(trace == min(trace(0.8e4:1.2e4))); %the minimum of the trace is the negative peak of the capacitative transient.

[pks,locs]=findpeaks(abs(diff(diff(trace))),'minPeakheight',7); %use the first peak of the second derivative to find onset of cap transient
intStart= locs(1);
hold on;
plot((intStart*Dx),trace(intStart),'r.','markersize',10);
intEnd=1.3e4; %manually selected.
steadyStateIndex=intEnd-500:intEnd+500;
steadyState=mean(trace(steadyStateIndex)); %value to be substracted for integration, only want to integrate the transient and not leak.
baseline= mean(trace(1:1e4));
plot(steadyStateIndex*Dx,trace(steadyStateIndex),'g');
xlim([8e3 intEnd+2000]*Dx);

area(((intStart:intEnd)*Dx),trace(intStart:intEnd),steadyState);

% Create a steady state substracted trace for integration;
SSsub=trace-steadyState;
integral=trapz((1:numel(intStart:intEnd))*Dx,(SSsub(intStart:intEnd)*1e-12));% Trapezoidal numerical integration of the capacitative transient after it has been steady state substracted.
capacitance=round(integral*1e12/DV,1);% capacitance in pF
Rm=round((DV/((steadyState-baseline)*1e-12))/(1e6),1); %Rm in megaohm.
hline(steadyState);
hline(baseline);

title([' Cm=' num2str(capacitance) ' pF, Rm= ' num2str(Rm) ' MOhm' ', Vstep=' num2str(DV*1e3) ' mV']);
obj=scalebar('XLen',0.01,'XUnit','s','YLen',20,'YUnit','pA')

clear fileNames dataMat
else 
end

%% Analyze spontaneous firing rate (on cell, whole cell)
fileNames=nameCellSort(spontLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
   % dataMat(:,i)=medfilt1(temp.y,10);
   
    
figure;plot(temp.y)
prompt = 'set threshold for pk detection (prominence, either negative or positive)         ';
thresh = input(prompt)
close gcf


if i==2
    figure;plot((1:numel(temp.y))*temp.dx*1e3, temp.y,'k')
xlim([0 2000])
obj=scalebar('XLen',100,'XUnit','ms','YLen',20,'YUnit','pA')
axis off
box off
end




% extract spike indices using peakfind
if thresh>0
[pks,locs]=findpeaks(temp.y,'minPeakprominence',thresh,'minpeakDistance',200);
else
[pks,locs]=findpeaks(temp.y*-1,'minPeakprominence',-thresh,'minpeakDistance',200);
end
% calculate mean firing rate and firing rate CV
if numel(locs)>1
% FR(i)=1/(temp.dx*(mean(diff(locs)))); %this only weorks if the cell is
% active a lot. otherwise it only takes the parts of the recording whre the
% cellis firing.

recDur=numel(temp.y)*temp.dx;
FR(i)=numel(locs)/recDur;
if numel(locs)>20
CV(i)=std(diff(locs))/mean(diff(locs));
end

else
    FR(i)=0;
    CV(i)=nan;
end


% calculate AP mean waveform and half width
spikeIndex=locs(5:end-5);
if numel(spikeIndex)>0
window=5*1e-3/temp.dx;% first term is window in ms
for j=1:numel(spikeIndex)
     mat4STAPre(j,:)= temp.y([(spikeIndex(j)-window):(spikeIndex(j)+window)]);
end

STATracePre(i,:)=nanmean(mat4STAPre,1);
timeSTA=(1:size(STATracePre,2))*temp.dx*1e3;
end
clear mat4STAPre




if i==2
    
prom= 5; % pA
[pks2,locs2]=findpeaks(medfilt1(temp.y,20),'minPeakprominence',prom,'minpeakDistance',50);

IPSCRate=numel(locs2)/(numel(temp.y)*temp.dx);
end
% 
figure
plot(medfilt1(temp.y,10))
hold on
plot(locs2,pks2+2,'ro')



end
%    
% CV
% FR
if exist('STATracePre')==1
if size(STATracePre,1)==3
%baseline=prctile(STATracePre(3,1:round(window)),5); % use the peak of the
%second derivative instead to find the inflection point.
secDer=diff(diff(STATracePre(3,:)));
[~, SDPkIdx]=max(secDer);
baseline=STATracePre(3,SDPkIdx);
peak=max(STATracePre(3,:));
figure;plot(timeSTA,STATracePre(3,:))
hline(peak,'--k')
hline(baseline,'--k')
halfMax=baseline + 0.5*(peak-baseline);



sub=STATracePre(3,:)-halfMax;
[~,idx1]=min(abs(sub(1:window)))
[~,idx2]=min(abs(sub(window:end)))
idx2=idx2+window;
idx1=idx1*temp.dx*1e3
idx2=idx2*temp.dx*1e3
width=idx2-idx1; % in ms
hold on
 plot([idx1 idx2],[halfMax halfMax],'k')
title({['Width at half max = ' num2str(idx2-idx1) ' ms'],['Firing rate = ' num2str(FR(3)) 'Hz, CV= ' num2str(CV(3))]})
xlabel('time (ms)')
ylabel('Vm (mV)')
else
end
else
    width=nan;
end

clear dataMat

%% Ih section
if sum(IhLogical)>0
    
    
fileNames=nameCellSort(IhLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    dataMat(:,i)=medfilt1(temp.y,10);
end

figure;plot(t,median(dataMat,2))
ylabel('current (pA)')
xlabel('time (s)')
clear dataMat

else
end

%% IF section  
% CHANGED TO ONLY ANALYZE THE FIRST 250 MS OF STEP.
fileNames=nameCellSort(IFLogical);

%first, extract the start of the step by looking at first derivative
%averages
for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    derivative(:,i)=diff(temp.y);
    figure;plot(temp.y)
end
meandiff=median(abs(derivative),2);
diffMAD=mad(meandiff(1:1e3));
diffBaseline=mean(meandiff(1:1e3));
threshCrossing=meandiff>(diffBaseline+(20*diffMAD));
StartStepIdx = find(threshCrossing, 1, 'first'); %find index of first crossing of 10 MADs of mean derivative of steps. This corresponds with the start of the step.
StartStepIdx=5e4; %override for when it's whacky and I can just set it manually
EndIdx=StartStepIdx+ (0.25/temp.dx)% Index of 250ms after step start. This is the region for counting spikes for IO curves regardless of step duration.



Dx=temp.dx;


for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    tempNotes=temp.WaveNotes;
    
    idxString=strfind(tempNotes,'a2:');
    semicolons=strfind(tempNotes,';'); % find idx of semicolons
    idxStringEnd= min(semicolons((semicolons-idxString)>0))-1; % find index of first semicolon after start of step value.
    
    Step(i)=str2num(tempNotes(idxString+3:idxStringEnd)); % 
    
    dataMat(:,i)=medfilt1(temp.y,10);
    
    %find spikes and count them.
    diffTemp=medfilt1(diff(temp.y),10);
    diffMADtemp=mad(diffTemp(1:1e3));
    [pks,locs]=findpeaks(diffTemp,'minPeakprominence',diffMADtemp*40,'minpeakDistance',200);
    
    numSpikes(i)=numel(locs(locs<EndIdx & locs>(StartStepIdx+100))); % only count first 250 ms
    
    FRSteps(i)=1/(temp.dx*(mean(diff(locs(locs<EndIdx)))));
    
end
FRSteps(isnan(FRSteps))=0; %get rid of the nans for plotting

figure;plot(Step,numSpikes,'linewidth',2)
hold on 
plot(Step,numSpikes,'k.','markersize',10)
ylabel('spike number')
xlabel('current injection (pA)')

% figure;plot(Step,FRSteps)
% hold on 
% plot(Step,FRSteps,'ko')
% ylabel('firing rate (Hz)')
% xlabel('current injection (pA)')

figure;plot(temp.y)
xlim([StartStepIdx-10e3 EndIdx+ 10e3])
hold on
plot(locs,pks,'ro')
vline(StartStepIdx,'--k')
vline(EndIdx,'--k')



%%simple 1D scatter plot
% 
% figure;
% hold on
% jitt=rand(numel(data),1)/3
% plot(1-1/6+jitt,data,'k.','markersize',30)
% plot([0.5 1.5], [mean(data) mean(data)],'k','linewidth',2)
% xlim([0 2])
% set(gcf,'color','white');
% set(gca,'FontSize', 18);
% box off
% xticks({})

CV
FR
capacitance
Rm

Step
numSpikes
IPSCRate

fclose all


